
*********************************
** Load correct visits effects:
	estimates use output/causal_visits
		global veffectp = _b[Np]
		global veffecte = _b[Ne]
*********************************


*****************************************
** Load relevant dataset

foreach type in income kids {

use "$datapath/intermediate/collect_pdouble_`type'.dta", clear

	gen id = 1 if raise=="Pe" 
	replace id = 2 if raise=="Pp"
	
	bys lno: gen Nlno = _N
	tab Nlno
	keep if Nlno==3
	distinct lno
	
*****************************************


*****************************************
** Keep libraries that solve correctly: (95.1% of libraries solve correctly)
	
	gen similar = (abs(Np-mNp)<.1 ) * (factor==10)
	egen SIMILAR = sum(similar), by(lno)
	tab SIMILAR
	keep if SIMILAR==1
	
	
*****************************************
** Additional outcome variables: librarian utility and visits

gen lutility= (thetap*mQp + thetae*mQe)/1000000
gen mQt = mQp + mQe

	** add in library visits effects
	
	tempfile cf
	save `cf'
	
	use $datapath/clean/main_data.dta, clear 
		keep if year == 2018
		keep lno visits
		
		merge 1:m lno using `cf'
		keep if _merge==3
		drop _merge
		gen double mvisits = visits + ((mNp-Np)*($veffectp) + (mNe-Ne)*$veffecte)
*****************************************

************************
** Drop libraries that solve to corners

	gen tag = mNp<1 | mNe<1
	egen TT = max(tag), by(lno)
	tab TT if raise=="Pe" & factor==10
	drop if TT==1
	

************************


*****************************************
** Keep averages and create ratios:

preserve 
	gen xc10 = CS if factor==10
	egen XC10 = mean (xc10), by(lno)	
	gen rCS = CS/XC10-1
	
	gsort lno factor raise
	gen rCS_ratio = rCS/rCS[_n-1] if raise=="Pp"

	collapse (mean) CS (p50) rCS_ratio, by(id factor raise)

		gen double x_CS = CS if factor==10
		egen double X_CS = mean(x_CS)
		gen r_CS = (CS-X_CS)/X_CS

		
	br raise factor r_CS
	
	keep if factor==20

		gen xCS = r_CS if raise=="Pe"
		egen XCS = mean(xCS)
		gen CS_ratio = r_CS / XCS if raise=="Pp"
		
		su rCS_ratio CS_ratio

restore	
**************************



****************************************************
** by income quintile

** Bootstrap for standard errors of CS ratios 

preserve 
	clear 
	set obs 1 
	gen v=. 
	save $datapath/temp/bs_CSratio.dta, replace 
restore 

forvalues k=1(1)100 {
preserve
	set seed `k'
	bsample, strata(ino) cluster(lno)
	
	tab ino 
	
	gsort lno factor raise
	
	collapse (mean) CS (semean) seCS = CS, by(factor raise ino)

	** CS at baseline prices
		gen double x_CS = CS if factor==10
		egen double X_CS = mean(x_CS), by(ino)
		gen r_CS = (CS-X_CS)/X_CS
	
	** Compare with raised prices
	keep if factor==20

		gen xCS = r_CS if raise=="Pe"
		egen XCS = mean(xCS), by(ino X)
		gen CS_ratio = r_CS / XCS if raise=="Pp"
		keep if raise=="Pp"

		keep CS_ratio ino
		
		twoway (connect CS_ratio ino) , scheme(lean2) 

		gen counter = `k'
		append using $datapath/temp/bs_CSratio.dta 
		sleep 300
		save $datapath/temp/bs_CSratio.dta, replace 
restore 	
}



preserve 
	use $datapath/temp/bs_CSratio.dta, clear 
	collapse (sd) sd = CS_rat, by(ino)
	tempfile sd 
	save `sd'
restore 



	
	
	collapse (mean) CS (semean) seCS = CS, by(factor raise ino)

		gen double x_CS = CS if factor==10
		egen double X_CS = mean(x_CS), by(ino)
		gen r_CS = (CS-X_CS)/X_CS

		
		
	keep if factor==20

		gen xCS = r_CS if raise=="Pe"
		egen XCS = mean(xCS), by(ino X)
		gen CS_ratio = r_CS / XCS if raise=="Pp"
		
		keep if raise=="Pp"
		keep raise CS_ratio ino
		
		merge m:1 ino using `sd'
		gen min = CS_rat  - 1.96*sd 
		gen max = CS_rat  + 1.96*sd 

		twoway (scatter CS_ratio ino) (rcap max min ino), scheme(lean2) xtitle(`type' quintiles) ytitle({&Delta}CS ratio) legend(off) ylabel(0(5) 25)
		graph export "$figpath/fig2_counterfactual_`type'.pdf", as(pdf) replace 
}


